clear all;

% parameters
alfa = 2/6;
d = 0.04;
theta = 0.75; 
R = 1.03; 
delta = 0.95; 
rho = delta * (1 - d - theta) + theta / R;  
a = (1-rho);
% steady state
qss = a / (1-rho) ; 

Hss = 1 / alfa / delta / d * (qss)^(alfa/(1-alfa));

kss = qss ^ (1/(1-alfa)); 

bss = theta * qss * Hss / R; 

hss = a*Hss;
% LTV loosening

theta_loosening = 0.9;

rho_loosening = delta * (1 - d - theta_loosening) + theta_loosening / R; 

qss_loosening = a / (1-rho_loosening);

Hss_loosening = 1 / alfa / delta / d * (qss_loosening)^(alfa/(1-alfa));
% LTV tightening

theta_tightening = 0.6;

rho_tightening = delta * (1 - d - theta_tightening) + theta_tightening / R; 

qss_tightening = a / (1-rho_tightening);

Hss_tightening = 1 / alfa / delta / d * (qss_tightening)^(alfa/(1-alfa));
% Learning parameters and initialization

%g = 0.033; % gain parameter 
g_high = 0.02; % gain parameter 

g_high_loose = 0.02; % gain parameter //top blue

g_high_tight = 0.02; % gain parameter //bottom blue


nperiod = 50; % no. of simulation 


hphigh = [
qss
qss
qss
qss
qss
qss
];

mhigh0 = qss; % 1; % 1.073;

np0 = size(hphigh,1); % length of HP growth data for initialization

m_high0 = zeros(np0, 1) ;

m_low0 = zeros(np0, 1) ;


for i = 1: np0
   
    if i == 1
        m_high0(i) = mhigh0 + g_high * (hphigh(i) - mhigh0) ; 

    else
        m_high0(i) = m_high0(i-1) + g_high * (hphigh(i)-m_high0(i-1)) ;         
    end
  
end

%%
% Learning dynamics
   
   % expectation change + loosening scenario

m_loosening_higher = zeros(nperiod,1);
hpgrowth_loosening_higher = zeros(nperiod, 1); 
m_loosening_higher(1) = m_high0(end)+0.018; 
H_loosening_higher = zeros(nperiod,1);

for i = 1: nperiod

    if i == 1

        hpgrowth_loosening_higher(i) = (1 - rho * m_high0(end-1))...
            / (1 - rho_loosening * m_loosening_higher(i));
        hp_loosening_higher(i) = qss_loosening * hpgrowth_loosening_higher(i);
        H_loosening_higher(i) = (1-d)*Hss_loosening+ hp_loosening_higher(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_loosening_higher(i) = a*H_loosening_higher(i);
        m_loosening_1(i) = m_loosening_higher(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_higher(i))^g_high_loose *(1-rho_loosening*1)^g_high_loose;
        m_loosening_2(i) = m_loosening_1(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_1(i))^g_high_loose *(1-rho_loosening*m_loosening_higher(i))^g_high_loose;
        m_loosening_3(i) = m_loosening_2(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_2(i))^g_high_loose *(1-rho_loosening*m_loosening_1(i))^g_high_loose;
        m_loosening_4(i) = m_loosening_3(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_3(i))^g_high_loose *(1-rho_loosening*m_loosening_2(i))^g_high_loose;

%    elseif i == 2

    else

                m_loosening_higher(i) =  m_loosening_higher(i-1) + ...
            g_high_loose * (hpgrowth_loosening_higher(i-1) - m_loosening_higher(i-1)) ;
        hpgrowth_loosening_higher(i) = (1 - rho_loosening * m_loosening_higher(i-1))...
            /(1 - rho_loosening * m_loosening_higher(i));
        hp_loosening_higher(i) = hp_loosening_higher(i-1) * hpgrowth_loosening_higher(i);
        H_loosening_higher(i) = (1-d)*H_loosening_higher(i-1)+ hp_loosening_higher(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_loosening_higher(i) = a*H_loosening_higher(i);
        m_loosening_1(i) = m_loosening_higher(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_higher(i))^g_high_loose *(1-rho_loosening*1)^g_high_loose;
        m_loosening_2(i) = m_loosening_1(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_1(i))^g_high_loose *(1-rho_loosening*m_loosening_higher(i))^g_high_loose;
        m_loosening_3(i) = m_loosening_2(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_2(i))^g_high_loose *(1-rho_loosening*m_loosening_1(i))^g_high_loose;
        m_loosening_4(i) = m_loosening_3(i)^(1-g_high_loose)/(1-rho_loosening*m_loosening_3(i))^g_high_loose *(1-rho_loosening*m_loosening_2(i))^g_high_loose;
    end
end


%%
   %  expectation change + tightening scenario

m_tightening_higher = zeros(nperiod,1);
hpgrowth_tightening_higher = zeros(nperiod, 1); 
m_tightening_higher(1) = m_high0(end)-0.011; 
H_tightening_higher = zeros(nperiod,1);

for i = 1: nperiod

    if i == 1

        hpgrowth_tightening_higher(i) = (1 - rho * m_high0(end-1))...
            / (1 - rho_tightening * m_tightening_higher(i));
        hp_tightening_higher(i) = qss_tightening * hpgrowth_tightening_higher(i);
        H_tightening_higher(i) = (1-d)*Hss_tightening+ hp_tightening_higher(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_tightening_higher(i) = a*H_tightening_higher(i);
        m_tightening_1(i) = m_tightening_higher(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_higher(i))^g_high_tight *(1-rho_tightening*1)^g_high_tight;
        m_tightening_2(i) = m_tightening_1(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_1(i))^g_high_tight *(1-rho_tightening*m_tightening_higher(i))^g_high_tight;
        m_tightening_3(i) = m_tightening_2(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_2(i))^g_high_tight *(1-rho_tightening*m_tightening_1(i))^g_high_tight;
        m_tightening_4(i) = m_tightening_3(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_3(i))^g_high_tight *(1-rho_tightening*m_tightening_2(i))^g_high_tight;

%    elseif i == 2

    else

                m_tightening_higher(i) =  m_tightening_higher(i-1) + ...
            g_high_tight * (hpgrowth_tightening_higher(i-1) - m_tightening_higher(i-1)) ;
        hpgrowth_tightening_higher(i) = (1 - rho_tightening * m_tightening_higher(i-1))...
            /(1 - rho_tightening * m_tightening_higher(i));
        hp_tightening_higher(i) = hp_tightening_higher(i-1) * hpgrowth_tightening_higher(i);
        H_tightening_higher(i) = (1-d)*H_tightening_higher(i-1)+ hp_tightening_higher(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_tightening_higher(i) = a*H_tightening_higher(i);
        m_tightening_1(i) = m_tightening_higher(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_higher(i))^g_high_tight *(1-rho_tightening*1)^g_high_tight;
        m_tightening_2(i) = m_tightening_1(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_1(i))^g_high_tight *(1-rho_tightening*m_tightening_higher(i))^g_high_tight;
        m_tightening_3(i) = m_tightening_2(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_2(i))^g_high_tight *(1-rho_tightening*m_tightening_1(i))^g_high_tight;
        m_tightening_4(i) = m_tightening_3(i)^(1-g_high_tight)/(1-rho_tightening*m_tightening_3(i))^g_high_tight *(1-rho_tightening*m_tightening_2(i))^g_high_tight;
    end
end


% counterfactual: policy change

% policy change: loosening; learning

m_constant_1 = zeros(nperiod,1);
hpgrowth_constant_1 = zeros(nperiod, 1); 
m_constant_1(1) = m_high0(end); 
H_constant_1 = zeros(nperiod,1);

for i = 1: nperiod

    if i == 1

        hpgrowth_constant_1(i) = (1 - rho * m_high0(end-1))...
            / (1 - rho_loosening * m_constant_1(i));
        hp_constant_1(i) = qss_loosening * hpgrowth_constant_1(i);
        H_constant_1(i) = (1-d)*Hss_loosening+ hp_constant_1(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_constant_1(i) = a*H_constant_1(i);
        m_constant_1_1(i) = m_constant_1(i)^(1-g_high)/(1-rho_loosening*m_constant_1(i))^g_high *(1-rho_loosening*1)^g_high;
        m_constant_1_2(i) = m_constant_1_1(i)^(1-g_high)/(1-rho_loosening*m_constant_1_1(i))^g_high *(1-rho_loosening*m_constant_1(i))^g_high;
        m_constant_1_3(i) = m_constant_1_2(i)^(1-g_high)/(1-rho_loosening*m_constant_1_2(i))^g_high *(1-rho_loosening*m_constant_1_1(i))^g_high;
        m_constant_1_4(i) = m_constant_1_3(i)^(1-g_high)/(1-rho_loosening*m_constant_1_3(i))^g_high *(1-rho_loosening*m_constant_1_2(i))^g_high;

%    elseif i == 2

    else

                m_constant_1(i) =  m_constant_1(i-1) + ...
            g_high * (hpgrowth_constant_1(i-1) - m_constant_1(i-1)) ;
        hpgrowth_constant_1(i) = (1 - rho_loosening * m_constant_1(i-1))...
            /(1 - rho_loosening * m_constant_1(i));
        hp_constant_1(i) = hp_constant_1(i-1) * hpgrowth_constant_1(i);
        H_constant_1(i) = (1-d)*H_constant_1(i-1)+ hp_constant_1(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_constant_1(i) = a*H_constant_1(i);
        m_constant_1_1(i) = m_constant_1(i)^(1-g_high)/(1-rho_loosening*m_constant_1(i))^g_high *(1-rho_loosening*1)^g_high;
        m_constant_1_2(i) = m_constant_1_1(i)^(1-g_high)/(1-rho_loosening*m_constant_1_1(i))^g_high *(1-rho_loosening*m_constant_1(i))^g_high;
        m_constant_1_3(i) = m_constant_1_2(i)^(1-g_high)/(1-rho_loosening*m_constant_1_2(i))^g_high *(1-rho_loosening*m_constant_1_1(i))^g_high;
        m_constant_1_4(i) = m_constant_1_3(i)^(1-g_high)/(1-rho_loosening*m_constant_1_3(i))^g_high *(1-rho_loosening*m_constant_1_2(i))^g_high;

    end
end




% counterfactual: policy change
% policy change: loosening; learning

m_constant_2 = zeros(nperiod,1);
hpgrowth_constant_2 = zeros(nperiod, 1); 
m_constant_2(1) = m_high0(end); 
H_constant_2 = zeros(nperiod,1);

for i = 1: nperiod

    if i == 1

        hpgrowth_constant_2(i) = (1 - rho * m_high0(end-1))...
            / (1 - rho_tightening * m_constant_2(i));
        hp_constant_2(i) = qss_tightening * hpgrowth_constant_2(i);
        H_constant_2(i) = (1-d)*Hss_tightening+ hp_constant_2(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_constant_2(i) = a*H_constant_2(i);
        m_constant_2_1(i) = m_constant_2(i)^(1-g_high)/(1-rho_tightening*m_constant_2(i))^g_high *(1-rho_tightening*1)^g_high;
        m_constant_2_2(i) = m_constant_2_1(i)^(1-g_high)/(1-rho_tightening*m_constant_2_1(i))^g_high *(1-rho_tightening*m_constant_2(i))^g_high;
        m_constant_2_3(i) = m_constant_2_2(i)^(1-g_high)/(1-rho_tightening*m_constant_2_2(i))^g_high *(1-rho_tightening*m_constant_2_1(i))^g_high;
        m_constant_2_4(i) = m_constant_2_3(i)^(1-g_high)/(1-rho_tightening*m_constant_2_3(i))^g_high *(1-rho_tightening*m_constant_2_2(i))^g_high;

%    elseif i == 2

    else

                m_constant_2(i) =  m_constant_2(i-1) + ...
            g_high * (hpgrowth_constant_2(i-1) - m_constant_2(i-1)) ;
        hpgrowth_constant_2(i) = (1 - rho_loosening * m_constant_2(i-1))...
            /(1 - rho_loosening * m_constant_2(i));
        hp_constant_2(i) = hp_constant_2(i-1) * hpgrowth_constant_2(i);
        H_constant_2(i) = (1-d)*H_constant_2(i-1)+ hp_constant_2(i)^(alfa/(1-alfa)) / (alfa*delta);
        h_consumption_constant_2(i) = a*H_constant_2(i);
        m_constant_2_1(i) = m_constant_2(i)^(1-g_high)/(1-rho_tightening*m_constant_2(i))^g_high *(1-rho_tightening*1)^g_high;
        m_constant_2_2(i) = m_constant_2_1(i)^(1-g_high)/(1-rho_tightening*m_constant_2_1(i))^g_high *(1-rho_tightening*m_constant_2(i))^g_high;
        m_constant_2_3(i) = m_constant_2_2(i)^(1-g_high)/(1-rho_tightening*m_constant_2_2(i))^g_high *(1-rho_tightening*m_constant_2_1(i))^g_high;
        m_constant_2_4(i) = m_constant_2_3(i)^(1-g_high)/(1-rho_tightening*m_constant_2_3(i))^g_high *(1-rho_tightening*m_constant_2_2(i))^g_high;

    end
end


font_size = 10;

Exp_loosening_higher = (m_loosening_higher-1) *100 ;
Exp_tightening_higher = (m_tightening_higher-1) *100 ;
Exp_constant_1 =  (m_constant_1-1)*100;
Exp_constant_2 =  (m_constant_2-1)*100;

Exp_loosening_higher_plot = [0 Exp_loosening_higher']';
Exp_tightening_higher_plot = [0 Exp_tightening_higher']';
Exp_constant_1_plot = [0 Exp_constant_1']';
Exp_constant_2_plot = [0 Exp_constant_2']';

Hconsumptio_loosening = (h_consumption_loosening_higher-hss)/hss*100;
Hconsumptio_tightening = (h_consumption_tightening_higher-hss)/hss*100;
Hconsumptio_constant_1 = (h_consumption_constant_1-hss)/hss*100;
Hconsumptio_constant_2 = (h_consumption_constant_2-hss)/hss*100;

Hconsumptio_loosening_plot = [0 Hconsumptio_loosening]';
Hconsumptio_tightening_plot = [0 Hconsumptio_tightening]';
Hconsumptio_1_plot = [0 Hconsumptio_constant_1]';
Hconsumptio_2_plot = [0 Hconsumptio_constant_2]';

for i = 1: nperiod
m_loosening_5YA(i) = ((m_loosening_higher(i) +m_loosening_1(i) +m_loosening_2(i) +m_loosening_3(i) +m_loosening_4(i))./5 -1)*100;
m_tightening_5YA(i)  = ((m_tightening_higher(i) +m_tightening_1(i) +m_tightening_2(i) +m_tightening_3(i) +m_tightening_4(i))./5-1)*100;
m_constant_1_5YA(i) = ((m_constant_1(i) + m_constant_1_1(i) + m_constant_1_2(i) + m_constant_1_3(i) + m_constant_1_4(i))./5 -1)*100;
m_constant_2_5YA(i) = ((m_constant_2(i) + m_constant_2_1(i) + m_constant_2_2(i) + m_constant_2_3(i) + m_constant_2_4(i))./5 -1)*100;

end

m_loosening_5YA_plot = [0 m_loosening_5YA]';
m_tightening_5YA_plot = [0 m_tightening_5YA]';
m_constant_1_5YA_plot = [0 m_constant_1_5YA]';
m_constant_2_5YA_plot = [0 m_constant_2_5YA]';


originhp = (1 - rho * m_high0(end-1))/ (1 - rho * m_tightening_higher(1));

HousePrice_loosening = (hp_loosening_higher-qss)/qss *1/2 *100;
HousePrice_tightening = (hp_tightening_higher-qss)/qss *100;
HousePrice_1 = (hp_constant_1-qss)/qss *100;
HousePrice_2 = (hp_constant_2-qss)/qss *100;

HousePrice_loosening_plot = [0 HousePrice_loosening]';
HousePrice_tightening_plot = [0 HousePrice_tightening]';
HousePrice_1_plot = [0 HousePrice_1]';
HousePrice_2_plot = [0 HousePrice_2]';

hpgrowth_loosening_higher_plot(1,1) = 0;
hpgrowth_tightening_higher_plot(1,1) = 0;
hpgrowth_constant_1_plot(1,1) = 0;
hpgrowth_constant_2_plot(1,1) = 0;


%back out 
HPrice_Loosen = HousePrice_loosening'/100 *qss +qss;
HPrice_Tighten = HousePrice_tightening'/100 *qss +qss;
HPrice_Const1 = HousePrice_1'/100 *qss +qss;
HPrice_Const2 = HousePrice_2'/100 *qss +qss;

for i = 2: length(HPrice_Loosen)-1
hpgrowth_loosening_higher_plot(i,1) = HPrice_Loosen(i,1)./HPrice_Loosen(i-1,1);
hpgrowth_tightening_higher_plot(i,1) = HPrice_Tighten(i,1)./HPrice_Tighten(i-1,1);
hpgrowth_constant_1_plot(i,1) = HPrice_Const1(i,1)./HPrice_Const1(i-1,1);
hpgrowth_constant_2_plot(i,1) = HPrice_Const2(i,1)./HPrice_Const2(i-1,1);
end

hpgrowth_loosening_higher_plot = (hpgrowth_loosening_higher_plot-1)*100;
hpgrowth_tightening_higher_plot = (hpgrowth_tightening_higher_plot-1)*100;
hpgrowth_constant_1_plot = (hpgrowth_constant_1_plot-1)*100;
hpgrowth_constant_2_plot = (hpgrowth_constant_2_plot-1)*100;

hpgrowth_loosening_higher_plot(1,1) = 0;
hpgrowth_tightening_higher_plot(1,1) = 0;
hpgrowth_constant_1_plot(1,1) = 0;
hpgrowth_constant_2_plot(1,1) = 0;

hpgrowth_loosening_higher_plot = [0 hpgrowth_loosening_higher_plot']';
hpgrowth_tightening_higher_plot = [0 hpgrowth_tightening_higher_plot']';
hpgrowth_constant_1_plot = [0 hpgrowth_constant_1_plot']';
hpgrowth_constant_2_plot = [0 hpgrowth_constant_2_plot']';





% Figure 5
figure 

plot(Exp_loosening_higher_plot(1 : nperiod+1),'-o','linewidth', 2)

hold on
plot(Exp_tightening_higher_plot(1:nperiod+1),'-^','linewidth', 2)
plot(Exp_constant_1_plot(1:nperiod+1),'--','linewidth', 2)
plot(Exp_constant_2_plot(1:nperiod+1),':','linewidth', 2)

title('1-year-ahead House Price Growth Expectation','fontsize', 14)
ylabel('Deviation from Original Steady State (p.p.)', 'fontsize', 14)

legend('LTV loosening', 'LTV tightening','LTV loosening (baseline)','LTV tightening (baseline)', 'fontsize', 9,'Location','Northeast')

set(gca, 'fontsize', 14); 
ylim([-3,5]);

xlim([0, nperiod]);

width=1040/2;
height=768/2;
dpi =70;
width_in = width / dpi;
height_in = height / dpi;
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [width_in, height_in]);
set(gcf, 'PaperPosition', [0, 0, width_in, height_in]);
set(gcf, 'PaperPositionMode', 'manual');

saveas(gcf,"LTV_HPExpectation_1Y.pdf",'pdf')



% Figure 4 right subplot
figure

plot(Hconsumptio_loosening_plot(1 : nperiod+1),'-o','linewidth', 2)

hold on
plot(Hconsumptio_tightening_plot(1:nperiod+1),'-^','linewidth', 2)
plot(Hconsumptio_1_plot(1:nperiod+1),'--','linewidth', 2)
plot(Hconsumptio_2_plot(1:nperiod+1),':','linewidth', 2)

title('Consumption of Housing Service','fontsize', 14)
ylabel('Deviation from Original Steady State (p.p.)', 'fontsize', 14)

legend('LTV loosening', 'LTV tightening','LTV loosening (baseline)','LTV tightening (baseline)', 'fontsize', 9,'Location','Northeast')

set(gca, 'fontsize', 14); 
ylim([-10,22]);

xlim([0, nperiod]);

width=1040/2;
height=768/2;
dpi =70;
width_in = width / dpi;
height_in = height / dpi;
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [width_in, height_in]);
set(gcf, 'PaperPosition', [0, 0, width_in, height_in]);
set(gcf, 'PaperPositionMode', 'manual');

saveas(gcf,"LTV_HousingConsumption.pdf",'pdf')


% Figure 4 left subplot
figure

plot(HousePrice_loosening_plot(1 : nperiod),'-o','linewidth', 2)

hold on
plot(HousePrice_tightening_plot(1 : nperiod),'-^','linewidth', 2)
plot(HousePrice_1_plot(1 : nperiod),'--','linewidth', 2)
plot(HousePrice_2_plot(1 : nperiod),':','linewidth', 2)

title('House Price','fontsize', 14)
ylabel('Deviation from Original Steady State (p.p.)', 'fontsize', 14)

legend('LTV loosening', 'LTV tightening','LTV loosening (baseline)','LTV tightening (baseline)', 'fontsize', 9,'Location','Northeast')

set(gca, 'fontsize', 14); 
ylim([-30,50]);

xlim([0, nperiod]);


width=1040/2;
height=768/2;
dpi =70;
width_in = width / dpi;
height_in = height / dpi;
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [width_in, height_in]);
set(gcf, 'PaperPosition', [0, 0, width_in, height_in]);
set(gcf, 'PaperPositionMode', 'manual');

saveas(gcf,"LTV_HPrice.pdf",'pdf')




